﻿// Copyright 2005, 2006 - Morten Nielsen (www.iter.dk)
//
// This file is part of SharpMap.
// SharpMap is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// SharpMap is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License
// along with SharpMap; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 

// SOURCECODE IS MODIFIED FROM ANOTHER WORK AND IS ORIGINALLY BASED ON GeoTools.NET:
/*
 *  Copyright (C) 2002 Urban Science Applications, Inc. 
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *  
 * updated by asida@centrum.cz for Petra.Gis
 * 
 *
 */
using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;
using Petra.Gis.ReferenceSystemServices;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{

  /// <summary>
  /// represents the Lambert Conformal Conic (lcc) one or two standard paralel projection.
  /// This version is enhandec by 3 extra params: cone kartographic pole lat, lon
  /// and gauss sphere standard parallel
  /// </summary>
  /// <remarks>
  /// <para>The Lambert Conformal Conic projection is a standard projection for presenting maps
  /// of land areas whose East-West extent is large compared with their North-South extent.
  /// This projection is "conformal" in the sense that lines of latitude and longitude, 
  /// which are perpendicular to one another on the earth's surface, are also perpendicular
  /// to one another in the projected domain.</para>
  /// </remarks>
  [DataContract]
  public class LambertConicConformalProjection : ProjectionTransform
  {
    double _falseEasting;   //projection false easting
    double _falseNorthing;  //projection false northing

    private double es;            //eccentricity squared
    private double e;             //eccentricity
    private double center_lon;    //longitude of projection center
    private double center_lat;    //latitude of projection center
    private double ns;            //ratio of angle between meridian
    private double f0;            //flattening of ellipsoid
    private double rh;            //height above ellipsoid
    private double k0;            //scale factor

    // equation variables
    private double _semiMajor;

    private double rh1;			// height above ellipsoid
    private double con;			// Sign variable
    private double ts;			// small t
    private double theta;		// angle
    private long flag;		  // error flag
    private double sinphi;  // sin value
    private double dX;      // difference to false easting
    private double dY;      // difference to false northing

    /// <summary>
    /// for oblique projection this transform must provide ellipsoid to gauss sphere transform
    /// </summary>
    private SphericConformalTransform _sct;


    #region Constructors



    /// <summary>
    /// xaml constr.
    /// </summary>
    public LambertConicConformalProjection()
      : base() { }




    /// <summary>
    /// Creates an instance of an LambertConformalConic2SPProjection projection object.
    /// </summary>
    /// <remarks>
    /// <para>The parameters must contain:</para>
    /// <list type="table">
    /// <listheader><term>Parameter</term><description>Description</description></listheader>
    /// <item><term>latitude_of_origin</term><description>The latitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
    /// <item><term>central_meridian</term><description>The longitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
    /// <item><term>standard_parallel_1</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is nearest the pole.  Scale is true along this parallel.</description></item>
    /// </list>
    /// <para>The optional parameters are:</para>
    /// <list type="table">
    /// <item><term>standard_parallel_2</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is furthest from the pole.  Scale is true along this parallel.</description></item>
    /// <item><term>false_easting</term><description>The easting value assigned to the false origin.</description></item>
    /// <item><term>false_northing</term><description>The northing value assigned to the false origin.</description></item>
    /// <item><term>scale_factor</term><description>Scale factor that reduces sphere size.</description></item>
    /// </list>
    /// </remarks>
    /// <param name="proj">ProjectedCoordinateSystem with mandatory projection parameters</param>
    public LambertConicConformalProjection(ProjectedCoordinateSystem proj)
      : this(proj, false) { }





    /// <summary>
    /// Creates an instance of an LambertConformalConic2SPProjection projection object.
    /// </summary>
    /// <remarks>
    /// <para>The parameters must contain:</para>
    /// <list type="table">
    /// <listheader><term>Parameter</term><description>Description</description></listheader>
    /// <item><term>latitude_of_origin</term><description>The latitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
    /// <item><term>central_meridian</term><description>The longitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined.</description></item>
    /// <item><term>standard_parallel_1</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is nearest the pole.  Scale is true along this parallel.</description></item>
    /// </list>
    /// <para>The optional parameters are:</para>
    /// <list type="table">
    /// <item><term>standard_parallel_2</term><description>For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is furthest from the pole.  Scale is true along this parallel.</description></item>
    /// <item><term>false_easting</term><description>The easting value assigned to the false origin.</description></item>
    /// <item><term>false_northing</term><description>The northing value assigned to the false origin.</description></item>
    /// <item><term>scale_factor</term><description>Scale factor that reduces sphere size.</description></item>
    /// </list>
    /// </remarks>
    /// <param name="proj">ProjectedCoordinateSystem with mandatory projection parameters</param>
    /// <param name="inverse">Indicates whether the projection forward (meters to degrees or degrees to meters).</param>
    public LambertConicConformalProjection(ProjectedCoordinateSystem proj, bool createInverse)
      : base(proj, createInverse){ }


    #endregion





    ///// <summary>
    ///// LambertConformalConicProjection axis list
    ///// </summary>
    //public override AxisList Axes
    //{
    //  get
    //  {
    //    //// this is default LambertConformalConic axis list with correct orientation
    //    //// this orientation must remain as defined because of projection equations
    //    //return new AxisList()
    //    //{
    //    //  new Axis("X", Axis.East),
    //    //  new Axis("Y", Axis.North)
    //    //};
    //    return _projCS.Axes;
    //    // todo: if projection is created from cache and
    //    // user touches this property then Axes must get from sourceCS or targetCS !!
    //  }
    //}





    protected override bool parametersValid()
    {
      return Parameters.ContainsKey("latitude_of_origin")
        && Parameters.ContainsKey("central_meridian")
        && Parameters.ContainsKey("standard_parallel_1");
    }




    protected override void updateTransform()
    {

      center_lat = Parameters["latitude_of_origin"] * base.AngularParameterUnit.RadiansPerUnit;
      center_lon = Parameters["central_meridian"] * base.AngularParameterUnit.RadiansPerUnit;

      double lat1 = Parameters["standard_parallel_1"] * base.AngularParameterUnit.RadiansPerUnit;
      double lat2;

      if (!Parameters.ContainsKey("standard_parallel_2"))
        lat2 = lat1;
      else
        lat2 = Parameters["standard_parallel_2"] * AngularParameterUnit.RadiansPerUnit;

      if (!Parameters.ContainsKey("false_easting"))
        _falseEasting = 0;
      else
        _falseEasting = Parameters["false_easting"] * DistanceParameterUnit.MetersPerUnit;

      if (!Parameters.ContainsKey("false_northing"))
        _falseNorthing = 0;
      else
        _falseNorthing = Parameters["false_northing"] * DistanceParameterUnit.MetersPerUnit;

      _semiMajor = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMajorAxis *
        this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      e = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.Eccentricity;
      es = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.e2;

      if (Parameters.ContainsKey("oblique_pole_latitude") && Parameters.ContainsKey("oblique_standard_parallel"))
      {
        //resolve oblique...
        if ((Parameters["oblique_pole_latitude"] * AngularParameterUnit.RadiansPerUnit) < ProjMath.HALF_PI)
        {
          _sct = SphericConformalTransform.Create(Projection.GeographicCoordinateSystem, AngularParameterUnit, IsInverse, Parameters);
          GeographicCoordinateSystem gaussSRS = (GeographicCoordinateSystem)_sct.Source;

          //use gauss sphere radius instead of ellipsoid
          _semiMajor = gaussSRS.Datum.Ellipsoid.SemiMajorAxis * gaussSRS.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

          e = gaussSRS.Datum.Ellipsoid.Eccentricity;
          es = gaussSRS.Datum.Ellipsoid.e2;
        }
      }

      //if scale factor is given no standard_parallel_2 must be applied !!
      //scale factor simulates 2 standard parallels
      if (Parameters.ContainsKey("scale_factor"))
        k0 = Parameters["scale_factor"];
      else
        k0 = 1;


      
      double sin_po;                  // sin value
      double cos_po;                  // cos value
      double con;                     // temporary variable
      double ms1;                     // small m 1
      double ms2;                     // small m 2
      double ts0;                     // small t 0
      double ts1;                     // small t 1
      double ts2;                     // small t 2

      /* Standard Parallels cannot be equal or on opposite sides of the equator
      ------------------------------------------------------------------------*/
      if (Math.Abs(lat1 + lat2) < ProjMath.EPSLN)
        throw new ArgumentException("Equal latitudes for St. Parallels on opposite sides of equator.");
            
      //center_lon = c_lon;
      //center_lat = c_lat;
      ProjMath.Sincos(lat1, out sin_po, out cos_po);
      con = sin_po;
      ms1 = ProjMath.Msfnz(e, sin_po, cos_po);
      ts1 = ProjMath.Tsfnz(e, lat1, sin_po);
      ProjMath.Sincos(lat2, out sin_po, out cos_po);
      ms2 = ProjMath.Msfnz(e, sin_po, cos_po);
      ts2 = ProjMath.Tsfnz(e, lat2, sin_po);
      sin_po = Math.Sin(center_lat);
      ts0 = ProjMath.Tsfnz(e, center_lat, sin_po);

      if (Math.Abs(lat1 - lat2) > ProjMath.EPSLN)
        ns = Math.Log(ms1 / ms2) / Math.Log(ts1 / ts2);
      else
        ns = con;

      f0 = ms1 / (ns * Math.Pow(ts1, ns));
      rh = this._semiMajor * f0 * Math.Pow(ts0, ns);
    }





    /// <summary>
    /// Converts coordinates from lon,lat to projection x,y.
    /// </summary>
    /// <param name="lonlat">coordinate in radians</param>
    /// <returns>coordinate in meters</returns>
    protected override void geographicToProjectedN(ICoordinate p)
    {

      // if oblique, then apply spheric conformal transform first
      if (_sct != null)
        _sct.Transform(p);

      //// get to radians
      //p.lon *= _geogUnit.RadiansPerUnit;
      //p.lat *= _geogUnit.RadiansPerUnit;

      con = Math.Abs(Math.Abs(p.lat) - ProjMath.HALF_PI);
      if (con > ProjMath.EPSLN)
      {
        sinphi = Math.Sin(p.lat);
        ts = ProjMath.Tsfnz(e, p.lat, sinphi);
        rh1 = _semiMajor * f0 * Math.Pow(ts, ns);
      }
      else
      {
        con = p.lat * ns;
        if (con <= 0)
          throw new ApplicationException();

        rh1 = 0;
      }

      theta = ns * ProjMath.AdjustLongitude(p.lon - center_lon);

      p.X = k0 * (rh1 * Math.Sin(theta) + _falseEasting);

      if(rh == 0)
        p.Y = k0 * (rh1 * Math.Cos(theta) + _falseNorthing);  //origin is at top of cone...
      else
        p.Y = k0 * (rh - rh1 * Math.Cos(theta) + _falseNorthing);

      //// convert to output units
      //p.X /= _projUnit.MetersPerUnit;
      //p.Y /= _projUnit.MetersPerUnit;


    }






    /// <summary>
    /// Converts x,y coordinate in projected units to lon,lat.
    /// </summary>
    /// <param name="p">Point in x,y</param>
    /// <returns>Transformed point as lon, lat</returns>
    protected override void projectedToGeographicN(ICoordinate p)
    {

      ////convert source coordinate to meters
      //p.X *= _projUnit.MetersPerUnit;
      //p.Y *= _projUnit.MetersPerUnit;

      p.X /= k0;
      p.Y /= k0;

      flag = 0;
      dX = p.X - this._falseEasting;

      if(rh == 0)
        dY = p.Y + this._falseNorthing; //origin is at top of cone...
      else
        dY = rh - p.Y + this._falseNorthing;

      if (ns > 0)
      {
        rh1 = ProjMath.Hypot(dX,dY);
        con = 1.0;
      }
      else
      {
        rh1 = -ProjMath.Hypot(dX, dY);
        con = -1.0;
      }
      theta = 0.0;
      if (rh1 != 0)
        theta = Math.Atan2((con * dX), (con * dY));

      if ((rh1 != 0) || (ns > 0.0))
      {
        con = 1.0 / ns;
        ts = Math.Pow((rh1 / (this._semiMajor * f0)), con);
        p.lat = ProjMath.Phi2z(e, ts, out flag);
        if (flag != 0)
          throw new ApplicationException();

      }
      else
        p.lat = -ProjMath.HALF_PI;

      p.lon = ProjMath.AdjustLongitude(theta / ns + center_lon);

      //// update to target units
      //p.lon /= _geogUnit.RadiansPerUnit;
      //p.lat /= _geogUnit.RadiansPerUnit;

      // if oblique
      if (_sct != null)
        _sct.Transform(p);

    }





    /// <summary>
    /// Returns the inverse of this projection.
    /// </summary>
    /// <returns>IMathTransform that is the reverse of the current projection.</returns>
    public override SRSTransform Inverse()
    {
      return new LambertConicConformalProjection(Projection, !IsInverse);
    }




    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }





    /// <summary>
    /// Oblique transform is sometimes called double transform and requres gauss sphere to put
    /// between projection and spheroid, so whole transform is composed as
    /// projection - gauss sphere - gauss sphere - reference ellipsoid.
    /// </summary>
    public override bool IsOblique
    {
      get 
      {
        if (_sct == null)
          return false;
        else
          return true;
      }
    }





    /// <summary>
    /// When LCC projection is oblique, this method turns it into non-oblique
    /// by developing LCC into 2 consecutive transforms (LCC and SphericConformal)
    /// </summary>
    /// <returns></returns>
    public override SRSTransform[] ToArray()
    {

      if (_sct == null)
        return base.ToArray();

      //a little trick,...clone original coordinate system and modify...
      ProjectedCoordinateSystem tmpPCS = (ProjectedCoordinateSystem)Projection.Clone();

      GeographicCoordinateSystem tmpGCS = tmpPCS.GeographicCoordinateSystem;

      //create geogr. coordinate system where gauss sphere is used as ellipsoid
      GeographicCoordinateSystem gauss = GeographicCoordinateSystem.CreateGaussSphereCS(
        tmpGCS, 
        Parameters["oblique_standard_parallel"] * AngularParameterUnit.RadiansPerUnit);

      //create spheric conformal transform as oblique
      //SphericConformalTransform tmpSC = new SphericConformalTransform(tmpGCS, tmpPCS.ProjectionDef.Parameters[], this.AngularParameterUnit, IsInverse);
      SphericConformalTransform tmpSC = SphericConformalTransform.Create(tmpGCS, AngularParameterUnit, IsInverse, tmpPCS.ProjectionDef.Parameters);
      
      //remove all oblique params from projection def
      tmpPCS.ProjectionDef.Parameters.Remove("oblique_standard_parallel");
      tmpPCS.ProjectionDef.Parameters.Remove("oblique_pole_latitude");
      tmpPCS.ProjectionDef.Parameters.Remove("oblique_pole_longitude");

      tmpPCS.GeographicCoordinateSystem = gauss;

      //and create new LCC from new definition as non oblique
      LambertConicConformalProjection tmpLCC = new LambertConicConformalProjection(tmpPCS, IsInverse);

      if(!this.IsInverse)
        return new SRSTransform[] { tmpLCC, tmpSC };
      else
        return new SRSTransform[] { tmpSC, tmpLCC };
    }





    public static string[] ProjectionNames
    {
      get { return new string[] { 
        "lcc", 
        "lambert conic conformal", 
        "lambert_conic_conformal", 
        "lambert_conic_conformal_1sp",
        "lambert_conic_conformal_2sp",
      }; }
    }







    /// <summary>
    /// creates new instance of LambertConicConformal projection transform from coordinate system definition.
    /// </summary>
    /// <param name="pcs">projected coordinate system that must be Lambert conic conformal</param>
    /// <param name="getAsInverse">return inverse transform = geog. to proj.</param>
    /// <returns>new instance of LambertConicConformalProjection class</returns>
    public static SRSTransform Create(ProjectedCoordinateSystem pcs, bool createInverse)
    {
      return new LambertConicConformalProjection(pcs, createInverse);
    }




  }


}

